Off— diagonal Ground State Properties of a ID Gas of Fermi Hard 

Rods 



F. Mazzanti 1 , G. E. Astrakharchik 2 , J. Boronat 2 , and J. Casulleras 2 
1 Departament de Fisica i Enginyeria Nuclear, Comte Urgell 187, 
Universitat Politecnica de Catalunya, E-08036 Barcelona, Spain and 
2 Departament de Fisica i Enginyeria Nuclear, Campus Nord B4-B5, 
Universitat Politecnica de Catalunya, E-08034 Barcelona, Spain 

Abstract 

A variational Monte Carlo calculation of the one-body density matrix and momentum distri- 
bution of a system of Fermi hard rods (HR) is presented and compared with the same quantities 
for its bosonic counterpart. The calculation is exact within statistical errors since we sample the 
exact ground state wave function, whose analytical expression is known. The numerical results are 
in good agreement with known asymptotic expansions valid for Luttinger liquids. We find that 
the difference between the absolute value of the bosonic and fermionic density matrices becomes 
marginally small as the density increases. In this same regime, the corresponding momentum 
distributions merge into a common profile that is independent of the statistics. Non-analytical 
contributions to the one-body density matrix are also discussed and found to be less relevant with 
increasing density. 
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I. INTRODUCTION 



Quantum one-dimensional (ID) systems of bosons and fermions have attracted great 
attention in the last years, both from the experimental and theoretical points of view p], \l , 
S 3, Is, 6]. The combined effect of quantum fluctuations and reduced dimensionality produces 
new and intriguing features different or not present in three-dimensional (3D) systems, as 
for instance the nonexistence of a true Bose-Einstein condensate in the homogeneous phase 
(not even at T = 0). 

As in the boson case, quasi-lD systems of fermions are experimentally realized confining 
to zero point oscillations the radial motion of a 3D trapped cloud of atoms in a largely 
elongated harmonic trap. This is done acting on the system with two perpendicular laser 
beams forming a two-dimensional optical lattice. In this setup, the transverse frequency 
ujj_ of the harmonic trap is much larger than the longitudinal frequency u\\, and such that 
the Fermi energy of the confined gas satisfies the condition Ep = Nhjj\\/2 <C hio^. Under 
these circumstances, an array of quasi-lD systems is created, thus allowing for a statistical 
evaluation of the most relevant quantities of interest p]]. 

Due to the Pauli exclusion principle, s-wave scattering between atoms of the same spin 
is not possible. In this way, the low-energy properties of the system are dominated by the 
interaction between atoms of different spin. In a pseudopotential description, the quasi-lD 
scattering length of the resulting interaction is directly related to the real 3D scattering 
length through the expression 

a 3 D \ a± 



a 1D = ^(l-C^\ , (1) 



with C = C(l/2)/v / 2 = 1.0326 and £(■) the Riemann zeta function 

A peculiarity of the one- dimensional world is that a certain number of exactly solvable 
many-body systems is known. The ground state energy of a system of bosons with 5- 
pseudopotential interactions (Lieb-Liniger model) has been obtained by using a Bethe 
ansatz solution in all interaction regimes, ranging from the weakly interacting Bose gas to 
the strongly correlated Tonks- Girardeau limit [lOj. The ground state energy of a system of 
particles interacting through a 1/z 2 potential (Calogero-Sutherland model) is also known 
for both fermions and bosons Contrary to the case of bosons, Fermi statistics preclude 
atoms from interacting via 5-pseudopotentials, which are commonly used to describe s-wave 
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scattering processes. Instead, ^-scattering is permitted between atoms of different spins in 
a system of two component fermions. Such systems have attracted great interest for a long 
time. There are two classes of excitations in homogeneous two-component ID Fermi gases, 
density and spin modes. In the case of attractive interactions sound waves propagate with 
a well defined velocity while spin waves exhibit a gap [ijj]. If the attraction is very strong, 
composite molecules consisting of two fermions of different spin are formed. 

Even if the energetic properties of a few ID systems can be determined exactly, much 
less is known about the correlation functions. Some short-distance properties (at least the 

n n n 

leading terms in a z — > expansion) are known for the one- |13[], two- |14{] and three- |15j 
particle correlation functions of the Lieb-Liniger model. The Tonks-Girardeau system is a 
special case where a finite number of terms in the expansion of the one-body density matrix is 
enough to provide a meaningful description of this quantity at all distances [3] . Moreover, it 



was noted in Refs. 
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171 ] that random matrix theory can be used to describe static density 



correlators in the Calogero- Sutherland model for certain values of the interaction parameter 
(A = 1/2,2 and 1, the later case being equivalent to the Tonks-Girardeau system). 

While the short-range properties of spatial correlators depend explicitly on the shape of 
the interaction potential, the long-range properties are quite generic and governed by the 
presence of phonons. The Luttinger model [3] describes universal long-range properties 
of all one- dimensional systems with a phononic (i.e. linear) excitation spectrum at low 
momenta. The long range description of the main ground state one- and two- particle 
correlation functions was obtained by Haldane in 1981 19|] . He described asymptotic series 
with universal power exponents (i.e. exponents that depend only on the density and the 
speed of sound) and non-universal series coefficients. All terms of the long-range expansion 
(with explicit expressions for the coefficients) of the one- 20] and two- j2l| particle correlation 
functions have been obtained using the replica method for the bosonic and the fermionic 
Calogero-Sutherland model, and for all possible interaction strengths including the Tonks- 
Girardeau limit. 

In spite of the progress achieved in analytical approaches, so far the only systematic way 
to obtain a complete description of the correlation functions relies on numerical methods. 



Recently, Monte Carlo methods have been success 



Lieb-Liniger |22J and the hard rod models 23j, |24| . Moreover, time-dependent Lieb-Liniger 



ully used to obtain these functions for the 



correlators have also been obtained using numerical summations of Bethe states 



231. In 
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much the same way, finite temperature correlators for the Lieb-Liniger model have been 



26]. 



calculated using density matrix renormalization techniques 

In this article we resort to Monte Carlo methods to analyze the most relevant ground 
state correlation functions of a single component system of fermionic hard rods. We note 
that while fermions with the same spin can not interact through a 5-pseudopotential due 
to the Pauli exclusion principle, a hard rod interaction is still permitted and provides the 
simplest interatomic potential defined by only one parameter a (the size of the hard-rod 
which equals its scattering length). 



II. RESULTS 



In a preceding paper 231 ] . the leading ground state properties of a Bose gas of ID hard 
rods were analyzed and discussed. In this work we extend this discussion to a population of 
fully polarized fermions interacting through the same spin-independent potential, namely 

+00 for \z\< a 
otherwise , 
corresponding to the many-body Hamiltonian 



V H r{z) 



(2) 
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2m / — ' dz 2 

j=i 3 



(3) 



i<3 



Despite the complexity of the interaction, both the ground state wave function and the 



energy are exactly known. The former is a Slater determinant of plane waves 
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=det 
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exp(ip k x k ) 



(4) 



2JV(JV-l)/2 
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vN\(L') N/2 f Kj 

where N is the total number of particles located in a box of length L with periodic boundary 
conditions. In this expression, L' = L — aN is the unexcluded length, while p' k = l-nn^jL 1 
with rtk an integer in the range —N to +iV plays the role of a single-quasiparticle momentum. 
Additionally, {x k = z k — (k — l)a} are a set of reduced coordinates for a given ordering of 
the true particle coordinates Zi < z 2 — a < z 3 — 2a < ■ ■ ■ < z N — a(N — 1). The ground state 
energy corresponding to this wave function reads 

E H r n 2 h 2 n 2 1 



Qm 



(5) 



with n = N/ L the linear density of the system. As happens in the 3D case of hard spheres, 
the scattering length of the hard rod potential equals the size of the rod, am = a. 

One important consequence of the constrains imposed by the restricted dimensionality 
in ID is the fermionization of bosons interacting through diverging potentials, where the 
strong repulsion between particles mimics the effect of the Pauli exclusion principle. In 
the particular case of hard rods, this duality is explicitly manifested in the form of \l/o, as 
in the Bose case the exact ground state wave function becomes the absolute value of the 
e„n in Eq. Q Q. I„ this case, therefore, a„ loca, quantltles depend exclusively o„ 
| \&o \ 2 , such as the energy, the static structure factor and the two-particle radial distribution 
function, are identical for bosons and fermions. However, other quantities not diagonal in 
configuration space are different due to the symmetry properties of \I/o- in this work, we 
analyze two of the most relevant non-diagonal ground state quantities, namely the one-body 
density matrix 



Both quantities have been evaluated via Monte Carlo sampling of the ground state wave 
function in Eq. (j3J) for a number of fermions N between 125 and 1001, located in a box of 
length L with periodic boundary conditions. Notice that since we use the analytical solution 
corresponding to the ground state of the Hamiltonian (j3J), the results of the simulation are 
exact in a statistical sense. 

The one-body density matrix satisfies the condition pi(0)/n = 1, which is a direct con- 
sequence of translational invariance and normalization. In the absence of a Bose-Einstein 
condensate, pi(z) decays to zero at large distances. For a system of fermions with a non- 
positive definite wave function, pi(z) changes sign with increasing distance. This is the case 
for instance of the ID free Fermi gas, where one has in the thermodynamic limit 



with an infinite number of nodes located at the points z m = m/n with m an integer. In the 
HR system at low densities, the average distance between particles is much larger than the 




(6) 



and its Fourier transform, the momentum distribution 




(7) 




(8) 
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rod size and the net effect induced by the potential is equivalent to a point-like boundary 
condition. In this limit \&o approaches the ground state wave function of the ID free Fermi 
gas (as can be easily checked from Eq. PJ setting a — > 0), and therefore p\{z) « pf^ G (z). 

The one-body density matrix p%(z) at the particle densities an = 0.1,0.2,0.4 and 0.6 
(pluses, crosses, stars and squares) is compared with pf^ G (z) (solid line) in Fig. dJ As 
it can be seen, little differences between pi(z) and PiE G { z ) arise at densities lower than 
an ~ 0.1. As the density increases, the main structure of the low-density pi(z) is kept while 
the strength at every point is depressed compared with p^ G {z). At intermediate and high 
densities na > 0.5 the oscillations are no longer visible at the scale of Fig. [TJ and essentially 
all of the strength of pi(z) is located around the origin. It is also worth to notice from the 
figure that the nodal structure of the one-body density matrix is poorly affected by particle 
correlations, keeping the nodes of pi(z) remarkably close to those of pf^ G (z). 

The one-body density matrix of any fermionic Luttinger liquid (as hard rods) admits the 
following asymptotic expansion valid at large distances z ^> n 19 1 

1 F m sin [27m (m + |) \z\] 

(n \z\) 1 / r ' 



Pl( Z ) 



77 



(9) 



(77 \z\Y m+ 2rv 

where 77 = 2K and K = nhn/mc is the Luttinger parameter written in terms of the sound 
velocity c. The coefficients F m and the value of 77 change with the interaction and are 
therefore model-dependent. For a free Fermi gas 77 = 2, F = 1/ir and F m>0 = 0. Apart 




anz 



FIG. 1: (Color online) p\{z) at the rod densities na = 0.1, 0.2, 0.4 and 0.6 (blue pluses, red crosses, 
green stars and blue squares) compared with the one-body density matrix of the ID free Fermi 
gas (solid line.) 
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from a few remarkable cases [la . l20j . these coefficients are, in general, unknown. For a 
system of hard rods, bosons or fermions, rj = 2(1 — an) 2 which is positive and decreases 
from 2 to with increasing density an G [0, 1]. High order terms in Eq. Q turn out to be 
then less relevant at large distances, and the asymptotic behavior of pi(z) when z — > oo is 
therefore dominated by the m = term. This term cancels at the nodes of p^£ G (z), that is, 
at the positions {z m = m/n}. At these points all other terms in the series vanish too, and 
the whole expression is zero. This fact explains why the nodes of pi(z) and pf£ G (z) are so 
close. Still, Eq. (jSJ) is an asymptotic expansion valid only beyond some healing distance, and 
therefore the first nodes of Pid{z) can deviate from those of PiE G { z )- This effect is small 
even at large densities and can hardly be appreciated in the an = 0.4 case of Fig. [TJ 

Another relevant aspect concerning the structure of the one-body density matrix of a 
fermionic system of hard rods is manifested when pi(z) is compared with Pib{z), the one- 
body density matrix of a system of boson hard rods of the same length and mass, and at the 
same density. Notice that Pib(z) is built as in Eq. ([6]) but starting from the corresponding 
bosonic hard rod ground state wave function, which is nothing but the absolute value of 



the wave function *&o(zi, z%, . . . , z N ) of Eq. (jl]) [23(. Figure [2] displays the absolute value of 
Pi(z) compared with Pib(z) (which is positive definite) for two particle densities an = 0.2 
and an = 0.6 (upper and lower panels, respectively). As it can be seen from the figure, both 
functions share a common short- distance behavior. It is easy to see from the symmetry of 
the one-body density matrix and the definition of the momentum distribution in Eq. (J7J) 
that the leading z — ► behavior of pi(z) and Pib( z ) is equal and given by 

?iM =1 _^l^ + cV), (10) 

n 6(1 — an) z 

where the coefficient of z 2 is proportional to the kinetic energy per particle which, for a 
system of hard rods (bosons or fermions), equals the total energy in Eq. (jSJ). The prediction 
of Eq. (1101) is shown as a dashed line in Fig. [2] for the two densities analyzed. Clearly, the 
quadratic approximation is too crude to describe p\ (z) but at very short distances. A better 
approximation is obtained when a truncated cumulant expansion is considered instead 



pi(*o r n2n2 ,2 



n 



exp 



z 



(11) 



6(1 — an) 2 

which is also displayed in the upper and lower panels of Fig. [2]with a solid line. As can be seen 
from the figure, the cumulant approximation works well both at low and high densities up to 




FIG. 2: (Color online) Absolute value of the one-body density matrix in logarithmic scale at the 
particle densities na = 0.2 and an = 0.6 (upper and lower panel, respectively). Black circles stand 
for fermions, red squares for bosons. The dashed and solid lines represent the approximations 
reported in Eqs. (I10p and (jlip . Inset: p\{z)/n compared with the truncated cumulant expansion 
of Eq. (|lip (solid and dashed lines, respectively). 

approximately the first node of pi(z). Beyond that point the approximation certainly breaks 
down because cumulant expansions can only be carried out on positive defined functions 
while pi(z) changes sign. In this sense, Eq. ( TTTT) is a better approximation to pis(z), although 
visible differences remain in the tails when plotted in logarithmic scale. At intermediate and 
high densities, however, there is no practical distinction between the fermionic and bosonic 
cases since the strength in pi(z) after the first oscillation is remarkably low. This is seen 
also in the inset of Fig. [2j where the dashed line is the cumulant approximation while p\ (z) 
and Pib(z) are displayed with a solid line and can not be distinguished from each other. 

At distances larger than the position of the first node, pi(z) and Pib{z) differ more 
significantly, specially at low densities. However, at large densities |pi(-z)| and Pib(z) are 
almost identical as seen from the figure, the main difference between the two functions being 
that P\b{z) is always positive while pi{z) changes sign each time a node is crossed. 



S 



This striking fact can be understood by direct inspection of the particle configurations 
contributing to the one-body density matrix. By definition, pi(z) is related to the probability 
of destroying a particle at the origin and creating a new one at a distance z, as given by 
(\Ert (2)^(0)) with \1/ and W field operators. In first quantization, p±(z) measures the overlap 
between the wave functions corresponding to different particle configurations where one 
particle has shifted its position by an amount z (see Eq. fl6])) while all other particles are 
kept in their original positions. Clearly, the different symmetry properties of the Fermi 
and Bose wave functions can make the contribution of these configurations to the one-body 
density matrix be of different sign. 

A remarkable property of ID systems like the one analyzed in this work is the fact that 
the different particle configurations can be classified in N\ disjoint subspaces, according to 
the ordering of the particles. In this way, one can attach labels 1 to AT to the particles of 
a given configuration, and build all other subspaces by sorting the particles in a different 
order. One can then move from the original subspace Qq to any other subspace Q a by 
transposition of particle coordinates. In the Bose case, the wave function is always positive. 
In the Fermi case, however, moving from one subspace to another implies a change in sign 
equal to (— l) p , where P is the parity of the permutation that leads to that subspace starting 
from the coordinate ordering of Qq. 

Now consider the high density limit for a system of hard rods. Given the wave function (j3J) , 
the most probable configuration is one where all particles are equally spaced. Each rod of 
size a has, in average, a space I = L/N available. At high densities, I is only slightly larger 
than a and there is not enough room for a particle to be placed between two other particles. 
In this sense, the most probable configuration will not contribute to Pi{z) when z is larger 
than I, as schematically represented in Fig. [3^. The only configurations in subspace flo 
contributing to pi(z) are those less probable were particles are more packed as in Fig. [3b, 
making room for the movement of the selected particle. Belonging to Qq, these configurations 
have always positive sign since no particle reordering is required, and contribute the same 
for bosons and fermions. Of course, this mechanism is only possible at short distances and 
that explains why the low z behavior of pi{z) and Pib{z) is equal. 

As the distance z increases, the contribution from Q to pi(z) becomes negligible and the 
main contributions come from movements involving a change from Q to different subspaces 
Q a . As before, and depending on the density, it may happen that the most probable con- 
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FIG. 3: (Color online) Three different configurations contributing to pi(z). (a) shows the most 
probable configuration where all particles are equally spaced. In (b) one particle is displaced a 
distance z, but there is not enough room for it to be placed between two other rods. The only 
configurations in f^o that allow this are ones were many particles have been shifted to close positions, 
and these are much less probable at high densities. In (c) the displacement z is large and that 
can only be done moving from f^o to another £l a , that is, changing the order of the particles with 
respect to the original ordering. 

figuration can not contribute due to the lack of available space, but there are configurations 
in other subspaces were only the movement of a few particles with respect to the equally 
spaced configuration is involved, as shown in Fig [3t. These are the most probable ones 
and produce the major contributions to the one-body density matrix. These contributions 
have a sign (— l) p for fermions and +1 for bosons. At high densities, only one subspace Q a 
contribute significantly, and therefore all contributions have the same sign. If this sign is 
negative, that leads to a total contribution to the one-body density matrix that is positive 
for bosons and negative for fermions but equal in absolute value. That explains the observed 
behavior shown in the lower panel of Fig. [2j 

It is also clear from the above arguments that this phenomenon does not happen at 
low densities, since in this case many different particle configurations have non-negligible 
contributions for a given z, and even at large distances configurations from different Q a 
subspaces contribute. In the Bose case all these configurations have the same sign, while 
in the Fermi one some are positive and others are negative, leading to a total contribution 
that is noticeably different for bosons and fermions. 

The form of the exact wave function is quite peculiar as the long-range decay of 
the two-body term sin(7nr/Z/) is very "weak". In this way, correlations between particles 
are "strong" even at large distances of the order of the box size. This situation should 
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be contrasted with two- or three- dimensional systems, where the two-body Jastrow factor 



is significantly different from 1 only at short distances 29]. As a consequence of strong 
correlations, density matrix p\{z) of the HR system does not saturate to a constant value, 
but decreases as z increases. This is the reason of the absence of a true BEC in the one- 
dimensional system of hard rods, even at zero temperature. Furthermore, it is also easy to 
understand that the decay is greatly enhanced in the high-density regime. Indeed, as we 
see, contributions to the off-diagonal \z\^> n -1 elements of the one-body density matrix at 
high density come mainly from the scenario described in Fig. 3c. The higher the density, the 
higher the cost of moving particles in order to make enough room for the displaced particle, 
and the smaller the weight of such a contribution. As a result p%(z) decays faster at higher 
densities. 

Being the one-body density matrix an even function of its argument, the momentum 
distribution (its Fourier transform, see Eq. ([71)) is also even n(— k) = n(k). Consequently 
and as originally pointed out by Sutherland 30], all odd moments of n{k) are zero. Moreover, 
(k 2n ) can be written in configuration space as the expectation value of (Y2j\d n ^ /dzj\ 2 ) once 
an integration by parts is carried out, and this is independent of the sign of the wave function. 
Consequently, all finite moments of n(k) are equal for Fermi and Bose hard rods. If pi{z) 
and Pib{z) were analytic functions of their arguments and the corresponding momentum 
distributions had non-divergent moments (k 2m ) < +oo for all integer m > 0, then pi(z) 
and p\B {z) would be equal, which is not the case. This can happen if only a finite number 
of moments of the momentum distributions exist, or if the one-body density matrices have 
non-analytic contributions. In the case of hard rods, both things happen. On one side 
the interatomic potential is a source of non-analyticity to the wave function, due to the 
excluded length which makes ^o(zi, . . . , zn) and all its derivatives be zero when two or 
more rods overlap, both for fermions and bosons. We thus see that for fixed coordinates z 2 
to zn, f Z2 ,-,z N {zi) = tya(zi,Z2, ■ ■ ■ ,zn) is a non-analytic function of zi, and one integrates 
fz 2 ,-,z N {z) to get the one-body density matrix. 

Notice the non-analyticity of pi(z) is an effect produced by the interatomic potential, and 

thus it is expected to apply both to the Bose and to the Fermi cases. In the an — > limit 

the potential reduces to a point-like interaction with well known properties. For bosons, 

□ 

the system enters the Tonks-Girardeau regime [10( • For fermions, the effect of the potential 
is already taken into account by the antisymmetry of the wave function, and the system 
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behaves as a ID free Fermi gas. In the later case, the corresponding momentum distribution 
is nf£ G (k) = 8(kF — k) and all moments exist. On the contrary, in the Tonks-Girardeau limit 
n(k) presents a l/\k\ 4 large-/c tail that makes all moments (k m ) with m > 3 diverge 13|. The 
same conclusions about the relevance of the statistics can be drawn when the ID free Bose 
gas, with n(k) oc 5(k) and all moments equal to 0, is compared with its fermionic counter par t 
(Fermionic Tonks Girardeau), which has a momentum distribution n(k) = l/[l + (fc/2) 2 ] 31]. 
Furthermore, the behavior of the momentum distribution of other exactly solvable ID models 
describing interacting Luttinger liquids can be related to the hard rod system though the 
value of the Luttinger parameter. For instance, in the Calogero-Sutherland model 11[ with 
an interatomic potential of the form 

A(A- iw 2 

V{z) = \ 2 ; , (12) 
L z sin {ttz/L) 

the momentum distribution for bosons and fermions decays to zero at large momentum 
according to a 1/ 1 k | 2+2A law [^(J, and thus only a finite number of moments of n(k) exist. 
Apparently, therefore, the inclusion of an interatomic potential in a ID system leads to a 
power-law decay in n(k) that makes it possible for the bosonic and the fermionic systems to 
have different momentum distributions while sharing all existing moments. Unfortunately, 
checking this feature from a Monte Carlo simulation is very difficult due to the enhancement 
of the statistical noise in the tail of n(k) when the order m of the moment (k m ) increases. 

In higher dimensions, the zero temperature momentum distribution n(k) of a normal 
Fermi liquid presents a gap at \k\ — kp that is inversely proportional to the effective mass 
at the Fermi surface 32]. This gap is also present in the ID Free Fermi gas. In ID, however, 
even weak particle correlations break down this picture. This is the case of Luttinger liq- 
uids 118 



where Haldane theory applies. The Fourier transform of the leading (m = 0) term 
in the asymptotic expansion reported in Eq. which oscillates with frequency kp = nn, 
reveals that n(k) is continuous at k = ±kp when I/77 + rj/A > 1. For hard rods, where 
rj = 2(1 — an) 2 and an e [0, 1], this condition is always met. In this way, the momentum 
distribution of the hard rod system of fermions is continuous at k = ±kp at all densities. 
Furthermore, the same analysis indicates that the slope of n(k) at these points is infinite 
and negative when I/77 + n/A < 2. Once again, for hard rods this condition implies the ex- 
istence of a threshold density n c a = 1 - \/l - V3 w 0.48 below which dn(k)/dk at \k\ = k F 
diverges. 
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FIG. 4: (Color online) Momentum distribution of the gas of Fermi hard rod at different densities 
(upper panel). From top to bottom at low k: an = 0.3 (black circles), an = 0.4 (red squares), 
an = 0.5 (green diamonds) and an = 0.6 (blue stars). The lower panel displays n(k) at na = 0.3 
for different number of particles: N = 125 (black diamonds), N = 251 (red squares), and N = 1001 
(blue circles). 

These two aspects are illustrated in the upper and lower panels of Fig. HI The upper 
plot shows the momentum distribution n(k) of the system at the particle densities an = 
0.3,0.4,0.5 and 0.6. As can be seen, the behavior of n(k) at \k\ = kp changes with the 
density, with steeper slope as the density decreases. These results are compatible with an 
infinite derivative at \k\ = kp for the densities an = 0.4 and an = 0.3. In any case, the 
present simulations have been carried out for an even number of particles in order to prevent 
the ground state from being degenerate. Since kp = -nn = nN/L, with a momentum spacing 
Ak = 2%/L imposed by the use of periodic boundary conditions, no point falls exactly at 
| A; | = kp. As a consequence, one can only investigate the behavior of the momentum 
distribution at the Fermi momentum by increasing the number of particles in the simulation. 
The lower plot in Fig. H] shows n(k) around kp for an = 0.3, as a function of the number of 



13 




k/k 

F 



FIG. 5: (Color online) Momentum distribution of the gas of Fermi (black squares) and Bose (red 
circles) hard rod at different densities. The blue dashed line in the lower panel shows the momentum 
distribution associated to the gaussian one-body density matrix of Eq. (fTT]) . 

particles in the simulation, for the three values N = 125,251 and 1001. As it can be seen 
from the figure, the points closer to kp corresponding to the highest N seem to confirm the 
theoretical prediction of a continuous n(k) at that point, while the results obtained at lower 
N, with a coarser spacing, could lead to the wrong conclusion that there is a gap in n{k) at 
k = ±kp. 

We conclude the results section comparing the momentum distribution of the Fermi and 
Bose systems of hard rods as a function of the density. Figure [5] shows n(k) at the rod 
densities an = 0.3,0.5 and 0.7 (upper, middle and lower panels, respectively). As expected, 
the momentum distributions are quite different at low densities, as in this limit the fermionic 
system approaches the ID free Fermi gas prediction while the bosonic one reproduces the 
Tonks-Girardeau limit. As the density increases, however, both functions are smeared out 
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and the differences become less relevant, to the point that both curves overlap and can 
not be distinguished from each other at the highest density considered an = 0.7. We have 
also checked numerically that the lowest order moments (k°) and (k 2 ), the former being a 
normalization factor and the later the kinetic energy per particle (which equals the total 
energy per particle for hard rods), are the same for bosons and fermions at the three densities 
reported in the figure. Additionally and for an = 0.7, Fig. [5] shows with a dotted line the 
Fourier transform of the gaussian one-body density matrix of Eq. (Hip, which is also gaussian. 

It is clear from Fig. [5] that the low density momentum distributions of both bosons 
and fermions are far from being gaussian, but also that the differences between n(k) and 
ns{k) reduce when the density increases. Furthermore, n(k) and ns{k) seem to approach 
a common profile that looks gaussian at low momenta, but this happens at even higher 
densities. This means that in the high density limit the contribution of the non-analytical 
parts of pi(x) is less relevant, and also that in this limit the coefficients in the Taylor 
expansion of the analytical part of p\{x) are mostly dominated by the kinetic energy per 
particle, according to the cumulant expansion of Eq. (TTTj) . It is also apparent from the figure 
that the differences between the statistics vanish, as a function of the density, before the 
common profile is approached. Still, one can not deduce from the simulation the analytical 
form of the large-/c behavior of the momentum distributions, and in particular whether a 
true gaussian is reached when an — > 1. 

III. SUMMARY AND CONCLUSIONS 

In summary, we have studied the one-body density matrix and momentum distribution 
of a Fermi gas of hard rods, comparing them to their bosonic counterparts. We find that 
Pi(z) has a nodal structure quite close to that of the ID free Fermi gas, in agreement with 
Haldane's theory of Luttinger liquids. We have also discussed the analytical properties of 
Pi(z) to find that non-analytical parts become less relevant as the density increases. Our 
numerical simulations confirm that the momentum distribution does not present a sharp 
gap at the Fermi surface, but has infinite derivative at k — ±kp when the density is lower 
than a critical value an c = 1 — a/2 — ~ 0.48. Furthermore, n{k) for fermions and 
bosons share the same lower order moments while showing an overall fairly different shape. 
These differences reduce when the density increases, and remarkably the Bose and Fermi 
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momentum distributions approach a common limit at high densities. 
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